clear all
% load productivities and F distributions (from COUNTERFACTUAL_A_gamma_0.m)

% Simulate (gamma = 1 x estimated value), general equilibrium, revenue
% neutral


global Dividend

Dividend=-642-37.5+7.3;
tau = 0.1;
load fitted_model
% Given Par2 optimal:
n = 50;
r = 0.036;
Par1 = zeros(29,1);
load ParEst_eh0_time0_Mex ParEst_eh0_time0_Mex
load ParEst_eh0_time1_Mex ParEst_eh0_time1_Mex
Par1(1:24) = ParEst_eh0_time0_Mex;
Par1(21) = 10^Par1(21);
Par1(22:24) = Par1(22:24)*1000;
Par1(25) = 500*ParEst_eh0_time1_Mex(1);
Par1(26) = ParEst_eh0_time1_Mex(2);
Par1(27) = ParEst_eh0_time1_Mex(3);
Par1(28) = ParEst_eh0_time1_Mex(4);
Par1(29) = ParEst_eh0_time1_Mex(5);
gamma =  Par1(25); a = Par1(22); b1 = Par1(23); b2 = Par1(24);
alphaf = Par1(17); betaf = Par1(18); alphai = Par1(19); betai = Par1(20); theta = Par1(21);
alphaf_2=Par1(26); betaf_2=Par1(27); alphai_2=Par1(28); betai_2=Par1(29);
load W_time1_age0_eh0_comp11.out
w_f_g=W_time1_age0_eh0_comp11; % Read wages, g and f distributions
wf_1=w_f_g(:,1);        % Wage in formal sector
wi_1=w_f_g(:,2);        % Wage in informal sector
wf_2=w_f_g(:,7);      % Wage in formal sector
wi_2=w_f_g(:,8);      % Wage in informal sector

wht_f_1=(wf_1(n)-wf_1(1))/(n-1);
wht_i_1=(wi_1(n)-wi_1(1))/(n-1);
wht_f_2=(wf_2(n)-wf_2(1))/(n-1);
wht_i_2=(wi_2(n)-wi_2(1))/(n-1);

p_f_1_est=p_f_1;
p_i_1_est=p_i_1;
p_f_2_est=p_f_2;
p_i_2_est=p_i_2;

% Generate log productivity main pctiles 
[pvector_f_1_est,Fvector_f_1_est]=main_quantiles2(1-Ff_1,p_f_1_est,n);
[pvector_i_1_est,Fvector_i_1_est]=main_quantiles2(1-Fi_1,p_i_1_est,n);
[pvector_f_2_est,Fvector_f_2_est]=main_quantiles2(1-Ff_2,p_f_2_est,n);
[pvector_i_2_est,Fvector_i_2_est]=main_quantiles2(1-Fi_2,p_i_2_est,n);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Solve for wage offer parameters

alphaf_old=alphaf;
betaf_old=betaf;
alphai_old=alphai;
betai_old=betai;

iter_wages=1;

for alphaf=alphaf_old*0.99%(0.98:0.01:1.02)
    for betaf=betaf_old*1%(0.98:0.01:1.02)
        for alphai=alphai_old*1%(0.98:0.01:1.02)
            for betai=betai_old*1.01%(0.98:0.01:1.02)
                
                Par1(17) = alphaf;
                Par1(18) = betaf;
                Par1(19) = alphai;
                Par1(20) = betai;
                
                [Ff_1,ff_1,Fi_1,fi_1,Ff_2,ff_2,Fi_2,fi_2]=solve_wage_offer_dist(Par1,wf_1,wi_1,wf_2,wi_2,n);
                [Wff,Wfi,Wfn,Wif,Wnf,Wii,Win,Wni,Wnn, Rw_fi_ni, Rw_ii_ni, Rw_ii_in, Rw_if_in]=solve_value_function_A3_C(Par1,Ff_1,wf_1,Fi_1,wi_1,Ff_2,wf_2,Fi_2,wi_2,n);
                [HHstates,wage_moments,transitions,gff, gfi, gif, gii, gfn, gnf, gin, gni]=data_simulation_A12_end(Par1,ff_1,wf_1,fi_1,wi_1,ff_2,wf_2,fi_2,wi_2,Wff,Wfi,Wfn,Wif,Wnf,Wii,Win,Wni,Wnn);
               
                % marginal distributions
                mff=HHstates(1);
                mfi=HHstates(2);
                mfn=HHstates(3);
                mif=HHstates(4);
                mnf=HHstates(5);
                mii=HHstates(6);
                min_=HHstates(7);
                mni=HHstates(8);
                mnn=HHstates(9);
                mf_1=mff+mfi+mfn;
                mi_1=mif+mii+min_;
                mf_2=mff+mif+mnf;
                mi_2=mfi+mii+mni;
                
                gf_1 = zeros (50,1); gf_2 = zeros (50,1); gi_1 =  zeros (50,1); gi_2 = zeros (50,1);
 for i = 1:50
 gf_1(i) =  (mff*(sum(gff(i,:))) + mfi*(sum(gfi(i,:))) + mfn*gfn(i))/(mff+mfi+mfn);   
 end

  for i = 1:50
 gi_1(i) =  (mif*(sum(gif(i,:))) + mii*(sum(gii(i,:))) + min_*gin(i))/(mif+mii+min_);   
  end
 
   for i = 1:50
 gf_2(i) =  (mff*(sum(gff(:,i))) + mif*(sum(gif(:,i))) + mnf*gnf(i))/(mff+mif+mnf);   
   end
 
    for i = 1:50
 gi_2(i) =  (mfi*(sum(gfi(:,i))) + mii*(sum(gii(:,i))) + mni*gni(i))/(mfi+mii + mni);   
    end
 
                    
                % firm size
                el_f_1=mf_1*gf_1./ff_1;
                % el_f_1(1)=el_f_1(2); % may need smoothing
                % el_f_1(n)=el_f_1(n-1);
                el_i_1=mi_1*gi_1./fi_1;
                % el_i_1(1)=el_i_1(2); % may need smoothing
                % el_i_1(n)=el_i_1(n-1);
                el_f_2=mf_2*gf_2./ff_2;
                % el_f_2(1)=el_f_2(2); % may need smoothing
                % el_f_2(n)=el_f_2(n-1);
                el_i_2=mi_2*gi_2./fi_2;
                % el_i_2(1)=el_i_2(2); % may need smoothing
                % el_i_2(n)=el_i_2(n-1);
                
                der_el_f_1=[Inf;el_f_1(2:n)-el_f_1(1:n-1)]/wht_f_1;
                der_el_i_1=[Inf;el_i_1(2:n)-el_i_1(1:n-1)]/wht_i_1;
                der_el_f_2=[Inf;el_f_2(2:n)-el_f_2(1:n-1)]/wht_f_2;
                der_el_i_2=[Inf;el_i_2(2:n)-el_i_2(1:n-1)]/wht_i_2;
                
                % productivity
                p_f_1=(wf_1+el_f_1./der_el_f_1)*(1+tau);
                p_i_1=wi_1+el_i_1./der_el_i_1;
                p_f_2=(wf_2+el_f_2./der_el_f_2)*(1+tau);
                p_i_2=wi_2+el_i_2./der_el_i_2;
                
                [pvector_f_1,Fvector_f_1]=main_quantiles2(1-Ff_1,p_f_1,n); 
                [pvector_i_1,Fvector_i_1]=main_quantiles2(1-Fi_1,p_i_1,n);
                [pvector_f_2,Fvector_f_2]=main_quantiles2(1-Ff_2,p_f_2,n);
                [pvector_i_2,Fvector_i_2]=main_quantiles2(1-Fi_2,p_i_2,n);
                
                crit_wages(iter_wages)=max([abs(pvector_f_1./pvector_f_1_est-1);abs(pvector_i_1./pvector_i_1_est-1);abs(pvector_f_2./pvector_f_2_est-1);abs(pvector_i_2./pvector_i_2_est-1)]);
                %crit_wages(iter_wages)=sum([(pvector_f_1-pvector_f_1_est).^2;(pvector_i_1-pvector_i_1_est).^2;(pvector_f_2-pvector_f_2_est).^2;(pvector_i_2-pvector_i_2_est).^2]);

                matrix_wages(iter_wages,:)=[alphaf,betaf,alphai,betai];
                iter_wages=iter_wages+1
                
                
            end
        end
    end
end

% Find solution and run for these values
iter_wages
min(crit_wages)
bestproduc=min(find(crit_wages<=min(crit_wages)))

alphaf=matrix_wages(bestproduc,1)
betaf=matrix_wages(bestproduc,2)
alphai=matrix_wages(bestproduc,3)
betai=matrix_wages(bestproduc,4)

Par1(17) = alphaf;
Par1(18) = betaf;
Par1(19) = alphai;
Par1(20) = betai;

[Ff_1,ff_1,Fi_1,fi_1,Ff_2,ff_2,Fi_2,fi_2]=solve_wage_offer_dist(Par1,wf_1,wi_1,wf_2,wi_2,n);
[Wff,Wfi,Wfn,Wif,Wnf,Wii,Win,Wni,Wnn, Rw_fi_ni, Rw_ii_ni, Rw_ii_in, Rw_if_in]=solve_value_function_A3_C(Par1,Ff_1,wf_1,Fi_1,wi_1,Ff_2,wf_2,Fi_2,wi_2,n);
[HHstates,wage_moments,transitions,gff, gfi, gif, gii, gfn, gnf, gin, gni]=data_simulation_A12_end(Par1,ff_1,wf_1,fi_1,wi_1,ff_2,wf_2,fi_2,wi_2,Wff,Wfi,Wfn,Wif,Wnf,Wii,Win,Wni,Wnn);
               
wage_moments_gamma_1_GE_rev_neut = wage_moments;
Rw_1_f_gamma_1_GE_rev_neut = Rw_fi_ni; Rw_1_i_gamma_1_GE_rev_neut = Rw_ii_ni;  Rw_2_i_gamma_1_GE_rev_neut = Rw_ii_in; Rw_2_f_gamma_1_GE_rev_neut = Rw_if_in;

save Rw_1_f_gamma_1_GE_rev_neut Rw_1_f_gamma_1_GE_rev_neut
save Rw_2_f_gamma_1_GE_rev_neut Rw_2_f_gamma_1_GE_rev_neut
save Rw_1_i_gamma_1_GE_rev_neut Rw_1_i_gamma_1_GE_rev_neut
save Rw_2_i_gamma_1_GE_rev_neut Rw_2_i_gamma_1_GE_rev_neut
HH_states_gamma_1_GE_rev_neut = HHstates;
save wage_moments_gamma_1_GE_rev_neut wage_moments_gamma_1_GE_rev_neut
save HH_states_gamma_1_GE_rev_neut HH_states_gamma_1_GE_rev_neut


mff=HHstates(1);
mfi=HHstates(2);
mfn=HHstates(3);
mif=HHstates(4);
mnf=HHstates(5);
mii=HHstates(6);
min_=HHstates(7);
mni=HHstates(8);
mnn=HHstates(9);
Un_Head = mnf + mni + mnn; Un_Spouse = mnn + min_ + mfn; Inf_rate = (mnn + min_ + mni + mii);
Agg_stocks_gamma_1_GE_rev_neut= [Inf_rate; Un_Head; Un_Spouse];
save Agg_stocks_gamma_1_GE_rev_neut Agg_stocks_gamma_1_GE_rev_neut
gf_1 = zeros (50,1); gf_2 = zeros (50,1); gi_1 =  zeros (50,1); gi_2 = zeros (50,1);
 for i = 1:50
 gf_1(i) =  (mff*(sum(gff(i,:))) + mfi*(sum(gfi(i,:))) + mfn*gfn(i))/(mff+mfi+mfn);   
 end

  for i = 1:50
 gi_1(i) =  (mif*(sum(gif(i,:))) + mii*(sum(gii(i,:))) + min_*gin(i))/(mif+mii+min_);   
  end
 
   for i = 1:50
 gf_2(i) =  (mff*(sum(gff(:,i))) + mif*(sum(gif(:,i))) + mnf*gnf(i))/(mff+mif+mnf);   
   end
 
    for i = 1:50
 gi_2(i) =  (mfi*(sum(gfi(:,i))) + mii*(sum(gii(:,i))) + mni*gni(i))/(mfi+mii + mni);   
    end

% have to construct the joint probabilities (for ex, gff(wf_1,wf_2)) or integrate over the
% marginal distributions gf_1 and gf_2, whichever is easier..

Welfare_total=r*(sum(sum(mff*Wff.*gff+mfi*Wfi.*gfi+mif*Wif.*gif+mii*Wii.*gii)))+r*sum(mfn*Wfn.*gfn+mnf*Wnf.*gnf+min_*Win.*gin+mni*Wni.*gni)+r*mnn*Wnn;
Welfare_formal1=(r*(sum(sum(mff*Wff.*gff+mfi*Wfi.*gfi)))+r*sum(mfn*Wfn.*gfn))/(mff+mfi+mfn);
Welfare_informal1=(r*(sum(sum(mif*Wif.*gif+mii*Wii.*gii)))+r*sum(min_*Win.*gin))/(mif+mii+min_);
Welfare_nonemployed1=(r*sum(mnf*Wnf.*gnf+mni*Wni.*gni)+r*mnn*Wnn)/(mnf+mni+mnn);
Welfare_formal2=(r*(sum(sum(mff*Wff.*gff+mif*Wif.*gif)))+r*sum(mnf*Wnf.*gnf))/(mnf+mif+mff);
Welfare_informal2=(r*(sum(sum(mfi*Wfi.*gfi+mii*Wii.*gii)))+r*sum(mni*Wni.*gni))/(mni+mii+mfi);
Welfare_nonemployed2=(r*sum(mfn*Wfn.*gfn+min_*Win.*gin)+r*mnn*Wnn)/(min_+mfn+mnn);
Welfare_gamma_1_GE_rev_neut=[Welfare_total;Welfare_formal1;Welfare_informal1;Welfare_nonemployed1;Welfare_formal2;Welfare_informal2;Welfare_nonemployed2]; % IN THOUSANDS, IF NORMALIZED
Budget_balance= -1306.51*Inf_rate-Dividend % the government pays 1306.51 per family/quarter for SP 
save Welfare_gamma_1_GE_rev_neut Welfare_gamma_1_GE_rev_neut

% Generate LOG WAGES by the F distribution 
[wvector_f_1_est,Fvector_f_1_w]=main_quantiles(1-Ff_1,wf_1,n);
[wvector_i_1_est,Fvector_i_1_w]=main_quantiles(1-Fi_1,wi_1,n);
[wvector_f_2_est,Fvector_f_2_w]=main_quantiles(1-Ff_2,wf_2,n);
[wvector_i_2_est,Fvector_i_2_w]=main_quantiles(1-Fi_2,wi_2,n);
mean_wf_1f=log(sum(wf_1.*ff_1));
mean_wi_1f=log(sum(wi_1.*fi_1));
mean_wf_2f=log(sum(wf_2.*ff_2));
mean_wi_2f=log(sum(wi_2.*fi_2));
wage_moments_under_F=[[wvector_f_1_est;mean_wf_1f],[wvector_i_1_est;mean_wi_1f],[wvector_f_2_est;mean_wf_2f],[wvector_i_2_est;mean_wi_2f]];
save wage_moments_under_f_gamma_1_GE_rev_neut wage_moments_under_F

% Transitions 
Trans_head = [transitions(3); transitions(4); transitions(1); transitions(5); transitions(2); transitions(6); transitions(13); transitions(14)];
Trans_spouse = [transitions(9); transitions(10); transitions(7); transitions(11); transitions(8); transitions(12); transitions(15); transitions(16)];
Transitions_gamma_1_GE_rev_neut = [Trans_head,Trans_spouse];
save Transitions_gamma_1_GE_rev_neut Transitions_gamma_1_GE_rev_neut

% Constructing produtivities for simulations

% marginal distributions
mf_1=mff+mfi+mfn;
mi_1=mif+mii+min_;
mf_2=mff+mif+mnf;
mi_2=mfi+mii+mni;

% firm size
el_f_1=mf_1*gf_1./ff_1; 
% el_f_1(1)=el_f_1(2); % may need smoothing
% el_f_1(n)=el_f_1(n-1);
el_i_1=mi_1*gi_1./fi_1;
% el_i_1(1)=el_i_1(2); % may need smoothing
% el_i_1(n)=el_i_1(n-1);
el_f_2=mf_2*gf_2./ff_2; 
% el_f_2(1)=el_f_2(2); % may need smoothing
% el_f_2(n)=el_f_2(n-1);
el_i_2=mi_2*gi_2./fi_2;
% el_i_2(1)=el_i_2(2); % may need smoothing
% el_i_2(n)=el_i_2(n-1);

% M (total employed workforce) in 2005 = 41441076
% Fraction of firms in the formal sector in 2002 includes self-emp. (ENAMIN/INEGI: 23.7%)
% Number of formal firms (IMSS registry, BOSCH AND CAMPOS-VASQUEZ 2005): 800,000
% Number of firms = 800000/0.237
M_N=41441076/(800000/0.237);

pi_f_1=(p_f_1_est-wf_1*(1+tau)).*el_f_1*M_N/0.237;
pi_i_1=(p_i_1_est-wi_1).*el_i_1*M_N/0.763;
pi_f_2=(p_f_2_est-wf_2*(1+tau)).*el_f_2*M_N/0.237;
pi_i_2=(p_i_2_est-wi_2).*el_i_2*M_N/0.763;

profitFormalFirms_1=1./(M_N*mf_1/0.237)*sum(pi_f_1(1:25).*ff_1(1:25)); % per worker 
profitInformalFirms_1=1./(M_N*mi_1/0.763)*sum(pi_i_1(1:18).*fi_1(1:18)); % per worker 
totalprofitFirms_1=mf_1*profitFormalFirms_1+mi_1*profitInformalFirms_1;  
profitFormalFirms_2=1./(M_N*mf_2/0.237)*sum(pi_f_2(1:25).*ff_2(1:25)); % per worker 
profitInformalFirms_2=1./(M_N*mi_2/0.763)*sum(pi_i_2(1:18).*fi_2(1:18)); % per worker 
totalprofitFirms_2=mf_2*profitFormalFirms_2+mi_2*profitInformalFirms_2;  
Profit_gamma_1_GE_rev_neut=[totalprofitFirms_1;profitFormalFirms_1;profitInformalFirms_1;NaN;totalprofitFirms_2;profitFormalFirms_2;profitInformalFirms_2]; % per worker
save Profit_gamma_1_GE_rev_neut Profit_gamma_1_GE_rev_neut

% normalized firm size
[elvector_f_1,Fvector_f_1]=main_quantiles2(1-Ff_1,el_f_1,n)
[elvector_i_1,Fvector_i_1]=main_quantiles2(1-Fi_1,el_i_1,n)
[elvector_f_2,Fvector_f_2]=main_quantiles2(1-Ff_2,el_f_2,n)
[elvector_i_2,Fvector_i_2]=main_quantiles2(1-Fi_2,el_i_2,n)
mean_el_f_1=sum(el_f_1.*ff_1);
mean_el_i_1=sum(el_i_1.*fi_1);
mean_el_f_2=sum(el_f_2.*ff_2);
mean_el_i_2=sum(el_i_2.*fi_2);
Firm_size_gamma_1_GE_rev_neut=[[elvector_f_1,elvector_i_1,elvector_f_2,elvector_i_2];[mean_el_f_1,mean_el_i_1,mean_el_f_2,mean_el_i_2]];
save Firm_size_gamma_1_GE_rev_neut Firm_size_gamma_1_GE_rev_neut
